Volume 112, Number 5, September-October 2007 

Journal of Research of the National Institute of Standards and Technology 



[J. Res. Natl Inst. Stand. Technol. 112,271-281 (2007)] 

Convective Instabilities in 
Two Liquid Layers 



Volume 112 



Number 5 



September-October 2007 



G, B. McFadden and 
S. R. Coriell 

National Institute of Standards 
and Technology, 
Gaithersburg, MD 20899-8910 

K. F. Gurski 

Department of Mathematics, 
George Washington Univesity 
Washington, DC 20052 

D. L. Cotrell 

Lawrence Livermore National 

Laboratory, 

Livermore, CA 94550 

mcfadden@nist.gov 
sam.coriell@nist.gov 
katliarine.gurski@nist.gov 
cotrell2@llnl.gov 



We perform linear stability calculations 
for horizontal fluid bilayers, taking into 
account both buoyancy effects and 
thermocapillary effects in the presence 
of a vertical temperature gradient. To 
help understand the mechanisms driving 
the instability, we have performed both 
long-wavelength and short-wavelength 
analyses. The mechanism for the large 
wavelength instabiUty is complicated, 
and the detailed form of the expansion is 
found to depend on the Crispation and 
Bond numbers. The system also allows a 
conventional Rayleigh-Taylor instability 
if heavier fluid overlies Hghter fluid, and 
the long-wavelength analysis describes 
this case as well. In addition to the 
asymptotic analyses for large and 
small wavelengths, we have performed 
numerical calculations using materials 
parameters for a benzene-water system. 
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1. Introduction 

The study of the stability of a fluid-fluid interface is 
important in a number of scientific and technological 
applications. In this paper we consider two immiscible 
fluid layers separated by a horizontal planar interface 
subject to a vertical temperature gradient. This problem 
has been well studied both theoretically and experimen- 
tally [1-6], and the effects of various driving forces on 
the stability of the system have been taken into account. 
Examples include the effects of buoyancy (natural or 
Rayleigh-Benard convection [7]), the effects of bulk 
density differences (Rayleigh-Taylor instabilities [8]), 
and the effects of surface tension gradients along the 
interface (Marangoni instabilities [1]), 



One of the earliest papers on the two-layer problem 
was by Zeren and Reynolds [9], who carried out 
numerical calculations, a long wavelength analysis, and 
experiments on the benzene-water system. Given the 
limited computational resources at that time, Zeren and 
Reynolds obtained reasonable numerical results. With 
the increased computational power that is now avail- 
able more accurate calculations may be performed, 
including the possibility of oscillatory (in time) modes 
that were ignored by Zeren and Reynolds, By using 
modem programs for symbol manipulation, the long 
wavelength analysis by Zeren and Reynolds can be 
extended to higher order as well. 

In this paper we examine the linear stability of hori- 
zontal fluid bilayers, taking into account both buoyancy 
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effects and thermocapillary effects in the presence of a 
vertical temperature gradient. The problem depends on 
a number of dimensionless parameters, including the 
ratio of material properties in each layer, the Rayleigh 
number (dimensionless temperature gradient), the 
Marangoni number (dimensionless force due to surface 
tension gradients), the Bond number (dimensionless 
buoyancy force at the interface), and the Crispation 
number (a measure of interface deformability). We find 
that in some situations the onset of instability can be 
due to a overstable mode that oscillates in time, with a 
finite critical wavenumber. In other circumstances the 
instability is stationary in time with a finite critical 
wavenumber. The stationary instabilities can extend to 
small wavenumbers, though the critical instability usu- 
ally does not occur at zero wavenumber. To help under- 
stand the mechanisms driving the instability we have 
performed both long-wavelength and short-wavelength 
analyses. The mechanism for the large wavelength 
instability is complicated, and the detailed form of the 
expansion is found to depend on the Crispation and 
Bond numbers. The system allows a conventional 
Rayleigh-Taylor instability if heavier fluid overlies 
lighter fluid, and this case is included in the results as 
well. In addition to the asymptotic analyses for large 
and small wavenumbers, we have performed numerical 
calculations using materials parameters for a benzene- 
water system. 

The paper is organized as follows. Governing equa- 
tions are given in the next section. The numerical pro- 
cedure is described briefly in Sec. 3. Linear stability 
results for various cases are presented next, including a 
comparison of the numerical results with large and 
small wavenumber expansions. A discussion is present- 
ed in Sec. 5, followed by conclusions in Sec. 6. An 
appendix contains a summary of the expansion results. 

2. Equations 

We consider a semi-infinite horizontal bilayer sys- 
tem, with vertical heating across the layers. The unper- 
turbed upper layer (denoted by a) extends over the 
interval 0<z<H„, and the unperturbed lower layer 
(denoted by p) extends over the interval -Hp<z<0. 
Without loss of generality we consider linear stability 
results for a two-dimensional system. The horizontal 
coordinate extends from - oo < x < ©o^ and the velocity 
u has components in the x and z directions denoted by 
u and TV, respectively. 



2.1 Governing Equations in the Bulk 

In each layer, we consider the Boussinesq equations 

Vu = 0, (1) 

p u, + p(u ■ V)u + Vp= fiV^u - pgz , (2) 

7; + (u- V)r=KV'r. (3) 

Here p is the pressure, fi is the dynamic viscosity, T 
is the temperature, k is the thermal diffusivity, g is the 
gravitational acceleration, / is the time, and z is the unit 
vector in the z-direction (anti-parallel to gravity). We 
assume ju and k are uniform in each layer, and also 
assume the density p is uniform in all terms except the 
gravitational term, where the density is given by 



p = p(i-77[r-rj). 



(4) 



in each layer. Here p is the density in each layer at the 
reference temperature Tj^, and the thermal expansion 
coefficient r/ is assumed to be uniform in each layer. 

2.2 Boundary Conditions 

The upper boundary atz=H„ and the lower bound- 
ary at z = - //^ are assumed to be isothermal with no- 
slip boundary conditions. The temperature is continu- 
ous across the interface. 



CT = 0, 



(5) 



where ITJ = T"-T^ denotes the temperature jump 
across the interface. The tangential velocity is assumed 
to satisfy the no-slip condition 



[u-tl = 0. 



(6) 



where t is any unit vector tangent to the interface. The 
stress boundary condition is 

IP u(u ■ ii - i>J = [T ■ i^l - yr A + Vs y, (7) 

where n is a unit normal vector to the interface, ^jk = 
-p5jk + ju(duj/dXk + dui,/dXj) is the stress tensor, y is 
the surface tension, AT is the curvature, v„ is the normal 
velocity of the interface, and V^ is the surface gradient. 
Here our sign convention is that the curvature AT is 
defined to be positive for a spherical inclusion of j8 
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phase. For example, in two dimensions with an inter- 
face z = h(x, /), the curvature is /C=-h^/[\ + h^f^, the 
interface velocity is v^ = hj[l + h^f^, and the surface 
gradient of the temperature-dependent surface energy 
7= '){T) is given by 



Vsr = r^J>mt, 



Vi+^ 



(8) 



where y^ = dyfdT and t is the unit tangent vector to the 
interface in the direction of increasing x. Here h^^ indi- 
cates the second derivative ofh, etc. 

The interface is a material surface, so that we have 






(9) 
(10) 



Under these conditions the left hand side of Eq. (7) 
vanishes. The continuity of heat flux gives 



dn 



0, 



(11) 



where k is the thermal conductivity and 3773^ = 
(7 - h^T^ ) / v^ + h^ is the normal derivative of the temp- 
erature field in each layer. 

2.3 Base State 

We linearize about a quiescent base state (also indi- 
cated by bars). The thermal field is 



T''iz) = T,+G„z, 



in the a layer, and 



T^iz) = T,+G,z, 



(12) 



(13) 



in the j3 layer, where T^ is the unperturbed interface 
temperature. The temperature gradients in the base state 
satisfy 



= kG-k,G, 



(14) 



The pressure field in the base state is hydrostatic, 
with 



dp'' _ dp^ 



dz 



dz 



(15) 



2.4 Dimensionless Parameters and Governing 
Equations 

Following the treatment by Zeren and Reynolds [9], 
we make the equations dimensionless based on a length 
scale given by the total depth d = Hf^-\- H^, a time 
scale based on the thermal time d^/K^, a velocity scale 
K^ Id, a temperature scale G^ d, and a pressure scale 
v^K^p^ Id^, These scales introduce the dimensionless 
parameters 



V =- 






1 






(16) 



_M« 



* ^a I * a r^* a * 

^ = — ' ^ ="^' ^ =7^' i" =—"' (17) 

K^ Kp (j^ jlp 



K„ v.Kg dy 



„ gppd' ^^ -rrGpd' 
Bo = — - — , Ma = — 



(18) 



f^^^^ 



and the geometrical parameter i=- H^IH ^. Here 
^a = Ma/Pa ^^^ ^^ = M^/p^ ^re the kinematic vis- 
cosities in each layer, Pr is the Prandtl number, Ra is the 
Rayleigh number, Cr is the Crispation number. Bo is 
the Bond number, and Ma is the Marangoni number. 
The minus sign in the definition of Ma is conventional, 
since for most materials, as for the benzene-water 
system, we have 77^ < 0. 

We assume a horizontal wavenumber a and a tempo- 
ral growth rate a= a,+ ia^ . The perturbed quantities 
(indicated by tildes) then satisfy 

zaw«+d)«=0, (19) 

Pr-^cjw« +iap''lp* =v*(w« -a'w«), (20) 

Pr-^cjd)"+A"/p*=v*(C-«'^")+^*Rar, (21) 

cjf" +G*d)" = K(r -a'r ), (22) 

for z > 0, and 

zaw^+d)f =0, (23) 



Pr ^ au^ + iap^ =^t~ ^^^^ ' 



(24) 



Pr-' cT<y'' + pf = ft^, - aW + Ra f , (25) 

(26) 



erf + ft)" =ff-a'f^ 



forz<0. 
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The boundary conditions at z = are 

f''+G'h = fl'+h, (27) 

u"-u'^=0, (28) 

ip"-p'')-BoCT-Xp*-\)h+a'CT-'h=2iiiar -a^), (29) 

ill'p.^-il^)+iail/d)"-d)'')-iaMaif" + G'h) =0, (30) 
af=oh, (31) 



(b'^ =ah, 

k'f" = ff . 



z z 



(32) 
(33) 



Critical conditions are often determined experimen- 
tally by varying the temperature gradient across the 
system. The temperature gradient G^ appears in the 
dimensionless parameters Ma and Ra. Zeren and 
Reynolds [9] introduce the parameter 



r = 



Ra _ -pglpd^ 



Ma 



(34) 



which is independent of the temperature gradient, and per- 
form calculations by varying Ma for a fixed value of P. 

3. Numerical Implementation 

We solve the eigenvalue problem that governs the 
linear stability of the system by using two complemen- 
tary procedures. In the first approach, the equations are 
discretized using pseudo-spectral Chebyshev colloca- 
tion, and the resulting generalized matrix eigenvalue 
problem is solved using the package RGG from the 
EISPACK software library [10]. For a discretization 
with A^ degrees of freedom, this routine produces 
approximations to the first A^ eigenvalues of the system. 
The second approach is to use the two-point boundary 
value solver BVSUP [11], coupled with the root fmder 
SNSQ [12], both from the SLATEC library [13], to 
implement a method described by Keller [14] to solve 
the eigenvalue problem. The BVSUP procedure pro- 
vides a very accurate solution for a given eigenmode 
provided a good enough initial estimate is available for 
the root-fmding procedure. The pseudospectral method 
is efficient for small values of A^, and is well-suited for 
searching parameter space to detect real and complex 
eigenvalues. Rather than performing fine grid calcula- 
tions with the pseudospectral procedure, however, the 
coarse grid results from the pseudospectral method are 



often used as initial guesses for the BVSUP code. 
Continuation from previous solutions also is used once 
an eigenmode has been identified. 

The BVSUP software works in a single domain, so 
we have mapped the two layers to a common domain 
by setting 



=1 ' 

\ -H„z/Hp 



forO<z<7f„, 
foT -H„ <z<0. 



(35) 



SO that 0<z<Hfy,in each layer. We then have 

J [ d/dz forO< z< //^, 

'd^~\ {\li)dldz for-H^ <z<0, (36) 

where l = -H^IH^ . To simplify the treatment of the 
problem, we also introduce an auxiliary ordinary differ- 
ential equation in z for the interface h , by setting 



dz 



(37) 



which allows us to avoid eliminating % as an unknown 
from the interface boundary conditions. 

4. Results 

In this section we present numerical results for the 
linear stability of the bilayer system. We provide a 
comparison with previous results of Zeren and 
Reynolds [9], who consider benzene overlying water. 
We use the thermophysical parameters given in their 
Table II for a temperature of 16 °C, and also consider 
their case dl = H^I{H^ + H^ = 0.4. Zeren and Reynolds 
consider both the case of heating from above (positive 
Ma), where the main driving force is the Marangoni 
effect and buoyancy effects are expected to be stabiliz- 
ing, and the case of heating from below (negative Ma), 
where both buoyancy effects and Marangoni effects 
can produce instabilities. 

If the system is heated from above, Zeren and 
Reynolds compute a positive critical Marangoni num- 
ber of 1486, with a critical wavenumber of a = 2.6; they 
consider only stationary modes with a; = 0. For a = 2.6 
we find Ma= 1468.4635, in fair agreement with their 
results. We find a critical wavenumber to three digits of 
a = 2.66, with Ma = 1466.8951. Ferm and Wollkind 
[15] also obtained similar agreement in a comparison 
with Zeren and Reynolds, although they chose different 
thermophysical properties for the benzene-water 
system. 
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Our computed neutral stability curves for this case 
are shown in Fig. 1. As discussed by Zeren and 
Reynolds, the critical mode is mainly driven by 
Marangoni convection. We show the two stationary 
modes that are given in the second figure of Zeren and 
Reynolds; the convection pattern in these modes has a 
single cell in each layer, with negligible surface defor- 
mation. The convection is more concentrated near the 
interface for the higher wavenumber mode. We also 
find an oscillatory mode (o;- i^ 0) at intermediate wave- 
numbers and Marangoni numbers. The oscillatory mode 
has a minimum near a = 5.5, with Ma = 22379.593 and 
<7i = ± 165.1. The surface deflection for this mode is 
negligible. The critical branch of the stationary mode 
does not extend to very small wavenumbers, but 
instead shows a large-Ma asymptote as the wave- 
number tends to the value a = 0. 12. 



10^ 







A/ 



10 



10 



10 



10 



Fig. 1. Marginal stability curves of Marangoni number Ma versus 
wavenumber a for a system of benzene over water with heating from 
above. Solid curves represent stationary modes with a,-= 0, and the 
dashed curve is an oscillatory mode with a^-^0. The circles denote 
the results of the large wavenumber expansion. 



If the system is heated from below, Zeren and 
Reynolds compute a negative critical Marangoni 
number of Ma = - 6068, with a critical wavenumber of 
a = 9. For a = 9.0 we find Ma = - 6014.4082. We show 
our computed marginal stability curves for this case 
in Fig. 2. We find that the critical mode is actually 
oscillatory, with Ma = -3146.7277 for a = 4.6 and 
o; = ± 64.93. The surface deflection for this mode is 
negligible. This mode exhibits a single convective cell 
in the upper layer, and two vertically-stacked cells in 
the lower layer of comparable strength. This mode 
merges with another stationary branch that has a long- 
wavelength asymptote; o;- tends to zero on the oscilla- 



tory branch as the modes merge. On the long-wavelength 
branch, for a = 0.001 we find Ma -- 19318.282. 
Analytic results for long-wavelength modes are 
described in the appendix; for this case the asymptotic 
result is Ma = - 19318.06, in excellent agreement with 
the numerical results. 




10 
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Fig. 2. Marginal stability curves for a system of benzene over water 
with heating from below. Solid curves represent stationary modes 
with a^-= 0, and the dashed curve is an oscillatory mode with a^-^ 0. 
The circles denote the limiting value of the small wavenumber 
expansion. 



There is a single convection cell in each layer for the 
long wavelength mode, with a significant interface 
deflection. Examination of the eigenmode shows that 
there is down-flow in the layer beneath the elevated 
portion of the perturbed interface. As discussed by 
Scriven and Stemling [16] for the single layer case, 
Marangoni modes and buoyant modes can be distin- 
guished by the direction of the vertical flow in the layer 
beneath an interface elevation: there is downflow 
beneath elevations for Marangoni flow, and upflow 
beneath elevations for buoyancy-driven flow. This is 
consistent with the Zeren and Reynold's interpretation 
of the long-wavelength mode as driven by Marangoni 
effects. The other two stationary modes in Fig. 2 with 
minima near a = 10 have negligible surface deflection. 
For both modes the lower layer is nearly isothermal 
with a weak unicellular flow, and the upper layer 
exhibits two vertically-stacked convective cells for the 
mode with Ma --6000, and three vertically- stacked 
convective cells for the mode with Ma --50,000. 
There is in fact an entire family of additional stationary 
buoyant modes, not shown in Fig. 2, with larger values 
of I Ma I, each containing multiple vertically- stacked 
cells in the two layers. 
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5, Discussion 

For the case of heating from below, as shown in 
Fig. 2 there is a small wavenumber instability that 
asymptotically approaches a finite Marangoni number 
as the wavenumber a tends to zero. To examine this 
behavior in more detail, we have performed additional 
calculations with F = 0, which eliminates the effects of 
buoyancy. Results are shown in Fig. 3 for both heating 
from below (Ma < 0) and for heating from above 
(Ma > 0). Suppressing buoyancy eliminates the oscilla- 
tory modes that prevail at intermediate wavenumbers for 
r = 0.142. For heating from below with our parameter 
values the stationary mode has a small-wavenumber 
limiting behavior that is given approximately by the 
expression 



-Ma: 



2.2691xlQ'(l+7.6Q26a') 

l-8.8850a' 



-O(a'), (38) 



(see Eq. (51) in the appendix). Note that in this approxi- 
mation there is a pole at a = 0.335. This is in approxi- 
mate agreement with the vertical asymptote that is 
obtained numerically for a = 0.268 for heating from 
either above or below. In fact, plotting instead 1/Ma 
versus a produces a single smooth curve crossing the 
ji;-axis with 1/Ma = at a wavenumber a = 0.268. 



mode. In Fig. 4 the upper curve corresponds to the data 
given in Table I. This curve asymptotically approaches 
a small wavenumber limit for Ma = -2.27 x 10"^. The 
lower curve corresponds to setting Bo = 0. The sym- 
bols on the curves correspond to numerical results, and 
the curves themselves correspond to analytical results 
from the small wavenumber approximation given in the 
appendix. The small wavenumber results depend 
strongly on both Bo and Cr. For the parameter values 
given in Table I, we have Cr = 2.1724 x 10"^ and 
Bo =1.1905, giving 



_Ma^ ^1^1::^^ = 2.2691 xlO\ 



c^Cr 



and for Cr = 2. 1724 x 10^ and Bo = 0, 



-Ma-^^ = 1.67190xl0'a'. 
c^Cr 



(39) 



(40) 



For Cr = and Bo = 1.1905, the resulting Marangoni 
number is positive (corresponding to heating from 
above), with 



Ma=^a-+£^ = ^:i^^^+583.1869, (41) 




10 



Fig. 3. Marginal stability curves for a system of benzene over water 
with buoyancy suppressed (F = 0). The circles denote the hmiting 
value of the small wavenumber expansion, and the squares denote 
the results of the large wavenumber expansion. 

To help further understand the low wavenumber 
instability that is observed numerically, we have also 
performed numerical computations and asymptotic 
expansions that illustrate the effects of the Bond 
number Bo and the Crispation number Cr for this 



and numerical results for this case are shown in Fig. 5. 



10% 



10^ 



^ 10^^ 
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^^^^ ^ / 
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-Cr>aBo>0 
-Cr>0. Bo-0 



Vo-^ 



10' 



10 



Fig. 4. Marginal stability curves for a system of benzene over water 
with buoyancy suppressed (r = 0). The circles indicate computation- 
al results obtained for the parameter values given in Table I, and the 
solid curve denotes the results of the corresponding small wave- 
number expansion. The diamonds indicate computational results 
obtained with Bo = 0, and the dashed curve denotes the results of the 
corresponding small wavenumber expansion. 
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Table I. Dimensionless parameters for the benzene (a layer) and 
water (jS layer) system at r= 16 °C, from Ref. [9] 



ratio of densities 


p' = Pa 


% 


0.886 


ratio of dynamic viscosities 


/ = Mg 


/M^ 


0.605 


ratio of thermal conductivities 


k* = k. 


fk^ 


0.274 


ratio of thermal diflfusivities 


^"=^a 


f^p 


0.730 


ratio of thermal expansion coefficients 


v" = rio 


JVi5 


7.06 


Bond number over Crispation number 


Bo/Cr 




5.48(10^) 


Inverse Bond number 


Bo-^ 




0.840 


Prandtl number 


Pr 




8 


Rayleigh number over Marangoni number 


r 




0.142 


dimensionless thickness ratio 


H^/Hp 




1.5 



10^ 



10^ 



^ 10 




ID- 



10' 
10"' 



Fig. 5. Marginal stability curves for a system of benzene over water 
with buoyancy suppressed (r = 0). The diamonds indicate computa- 
tional results obtained with Cr = 0, and the solid curve denotes the 
results of the corresponding small wavenumber expansion. 

Here the coefficients Ci, C2, C4, and c^ in these expres- 
sions generally depend on the layer geometry and the 
remaining material constants and are given in the 
appendix; here we have evaluated these expressions for 
the values given in Table I. 

5.1 Large Wavenumbers 

The temperature and flow field for the Marangoni 
instability can be studied in the large wavenumber limit 
that is summarized in the appendix. We ignore buoyan- 
cy by setting F = 0, and also ignore interface deforma- 
tion by setting Cr = 0, which the numerical results indi- 
cate are good approximations in this case. In the large 
wavenumber limit the flow is confined to the vicinity of 
the interface and the effect of the upper and lower 
boundaries ^tz = H„ and z = - //^ is negligible. To help 
visualize the flow we introduce a two-dimensional 
streamfunction y/ with w = y/^ and u = -y/^ . Contours 
of the temperature and streamfunction near the inter- 



face are shown in Fig. 6. Here we have exaggerated the 
size of the perturbation to emphasize the distortion 
of the isotherms near the interface. To make the plots 
easier to interpret we have assumed equal thermal con- 
ductivities (k* = I), which equalizes the unperturbed 
temperature gradients in the two layers and facilitates 
comparison of the perturbed temperature fields in each 
layer. There is a temperature gradient along the inter- 
face (atz = 0), and the streamlines of the flow are along 
the interface. Over the single period of the flow shown 
in Fig. 6, there are four convective cells near the inter- 
face that alternately compress and expand the isotherms 
near the interface. For our parameters with k^ < 1, the 
Marangoni instability occurs for heating from above 
(Ma > 0) in the large wavenumber limit. If k^ > 1 the 
instabilities occur for heating from below (Ma < 0) 
instead. 




Fig. 6. Streamfunction contours (light lines) and temperature con- 
tours near the interface for the large wavenumber Marangoni insta- 
bility with a=l, k*=\.0. and k** = 0.5. The magnitude of the 
perturbation is exaggerated to emphasize the deformation of the 
temperature contours. 



5.2 Rayleigh-Taylor Instability 

The two-layer system can exhibit the classical 
Rayleigh-Taylor instability if heavier fluid overlies 
lighter fluid. In the simplest case the Rayleigh-Taylor 
instability can be understood by a simple potential 
energy argument that balances the increased surface 
energy of a deformed interface y = h{x) against the 
change in the gravitational potential energy of the 
displaced fluid, 



SiPa-p^)h=-7h^. 



(42) 
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In terms of our dimensionless variables, this takes the 
form 



-Bo=- 



(1-p*) 



(43) 



which can be seen as a factor in the normal stress 
balance boundary condition (29). 

In the situation we have studied above, we have 
lighter benzene overlying heavier water, so the 
Rayleigh-Taylor instability does not occur. To study 
this instability for our system with a minimal change in 
notation, we temporarily choose to change the direction 
of gravity while keeping benzene and water in the 
original orientation, so that the water and benzene 
are unstably stratified with respect to gravity. We take 
G^ < 0, so that buoyancy has a stabilizing effect on the 
system; the resulting sign conventions produce Ma < 0, 
Ra > 0, r < 0, and Bo < 0. In Fig. 7 we show the corre- 
sponding numerical results for the Rayleigh-Taylor 
instability. The dotted curve in Fig. 7 shows the curve 
-Bo = aV(l - p*) that holds for Ra = Ma = 0. The solid 
curve shows numerical results for Ra = 1.0 x 10^ and 
Ma = 0. The stabilizing effect of buoyancy is evident at 
small wavenumbers, where the system is then stable if 
I Bo I is sufficiently small. The marginal stability curve 
asymptotes at small wavenumbers to the value 



-c.CrRa 
Bo = -^ -6.45 

q(l-p*) 



(44) 



in our case, which follows from Eq. (45) in the appen- 
dix. The dashed curve shows numerical results for 
Ra = and Ma = - 1.0 x 10^ The destabilizing 
Marangoni effect in this configuration is evident at 
small wavenumbers. and there is a cut-off wavenumber 
at a = 0.333 where Bo = 0; at lower wavenumbers this 
mode merges into the Marangoni instability with posi- 
tive Bond numbers. This value of the cut-off wavenum- 
ber is too large to be quantitatively approximated by 
the small wavenumber expansion given by Eq. (51) in 
the appendix. For the benzene-water system with 
Ra/Ma = - 0.142, numerical results are qualitatively 
similar to those for the case with Ra = 0. 

Our small wavenumber expansion of the instability, 
as summarized in the appendix, corrects and extends a 
similar analysis given by Zeren and Reynolds [9], who 
note that at leading order in a only the interface posi- 
tion and temperature field are perturbed, with no lead- 
ing order velocity profiles. However, the asymptotic 
values for the critical Marangoni number do depend on 
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Fig. 7. Marginal stability curves for a water-benzene system that 
is unstably stratified with respect to gravity The dotted curve 
(Ma = Ra = 0) corresponds to the classical Rayleigh-Taylor instabil- 
ity in the absence of buoyancy given by- Bo = a /(I - p*). The solid 
curve represents numerical resuhs that include the effects of buoy- 
ancy with Ra= 1.0 X 10^ and Ma= 0, and the circles are the small 
wavenumber asymptotic results for this value of Ra. The dashed 
curve represents numerical results that include the Marangoni effect, 
with Ra = and Ma = - 1.0 X 10^ 

the viscosity of the fluid, so flow does have an effect on 
the stability of the system, even at small wavenumbers. 
The small wavenumber analysis proceeds by an expan- 
sion in powers of a^, and the critical Marangoni number 
is determined at a later point in the expansion where the 
0(a^) flow perturbation is significant. We also note that 
the expansion procedure is based on a limit of small but 
non-zero wavenumbers. In the special one-dimensional 
case a = the perturbed interface is flat and mass con- 
servation restricts its location to z = 0. In a sense this 
represents a discontinuity in the problem in the limit 
a ^ 0, since in the expansion procedure a shift in the 
interface location is obtained at leading order. This is 
somewhat analogous to the dependence of a perturbed 
interface of the form z=Hq cos a x, which conserves 
mass for a^O (i.e., the average interface position is 
zero), but represents a shifted flat interface at z = T/q for 
a = 0. In general, for real experiments the smallest 
allowable value of a will be limited by the finite later- 
al extent of the system, but whether or not a perturba- 
tion with a = is physically relevant depends on the 
specific problem under consideration. 

6. Conclusions 

We have performed linear stability calculations for 
horizontal fluid bilayers, taking into account both 
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buoyancy effects and thermocapillary effects. We con- 
sider the case of a lighter layer overlying a heavier 
layer, so that the base state is stably stratified in this 
sense. We find that the system can be linearly unstable to 
either heating from above (Ma > 0) or below (Ma < 0). 
The mechanism for the large wavelength (small a) 
instability is complicated (see Figs. 4 and 5). The 
detailed form of the expansion, viz, the exponent n in 
the leading order expansion Ma ~ a"" depends on the 
Crispation and Bond numbers. 

7. Appendix 

Here we consider the limits of large and small wave- 
numbers. We introduce the dimensionless layer widths 
H^ = H^/d= 1/(1 - i) and Hp = H^ld = - 1/(1 - I), 
where d = H^ + H^ and l=-H^ IH^ = - 2/3 in our 
calculations. 

7.1 Small Wavenumbers 

By performing long-wavelength asymptotics for 
Bo 9^ and Cr ;^ we find the leading order result 



Ma = 



c^Boip -1) 



(45) 



For the values given in Table I, we have 
Ci =11.4717, c^ =31.5840, c, =38.8300. (49) 

We note that in general q is positive, but C2 can 
change sign depending on the magnitude of ju% H^ , and 
H^=\-H^, For a fixed layer geometry (given H^ ), 
large and small values of//* thus produce a sign change 
in the asymptotic value of Ma for small a and F = 0; 
this corresponds to the dominant flow occurring prima- 
rily in either the upper or lower layer as the correspon- 
ding viscosity contrast is large. The expression (48) for 
C5 is more complicated and also can have either sign, 
depending on the parameter values. 

7.1.1 Case for Ra = 

If we set Ra = the linear eigenfunctions for the 
general problem can be written down explicitly, and the 
dispersion relation for (7= can be evaluated symboli- 
cally in closed form, though the resulting expression is 
lengthy and difficult to interpret. However, the result 
can be expanded for small wavenumbers, and produces 
expressions that are not too unwieldy. 

The general expression for (7=0 and Ra = can be 
found by expanding in the wavenumber a to give 



where 



= fl'°{c,[p*-l]Bo-C2CrMa} 



(50) 



q =80#„#,[^„ +k'H^MH^ +M^^,], (46) W^{-Ma(c3Ci^cJpM]Bo)-(q -cJpM]Bo)} +0(a'^) 



c,=m[Hl-^'Hl] 



c.^^^^iUHlHlVn ^* pH^-k'H^] 



(47) 



(48) 



+3[npHl -k'^'H'A+UH^HJiiplf^ -kym 



The expression can be solved for Ma to give a ration- 
al expression of the form 



Ma 



q[p*-l]Bo-(q-c,[p*-l]Bo)a^ 
C2Cr+(c3Cr+C4[p* -l]Bo)a' 



+0{a\ (51) 



The constants Ci and C2 are as given above, and the 
additional constants C3, C4, and c^ are given by 



This expression is in good agreement with the 
numerical results for small wavenumbers as shown in 
Fig. 2. For Ra = 0, this expression agrees with that 
given by Zeren and Reynolds [9] up to a difference in 
sign; their expression also agrees with their numerical 
results for Ra = if the sign of their expression is 
changed. For Ra ;^ our expression gives results which 
differ in both magnitude and sign from that given by 
Zeren and Reynolds. We note that the leading order 
result in Eq. (45) is proportional to Bo and inversely 
proportional to Cr, so that the special cases Bo = or 
Cr = require additional attention; we consider these 
cases in more detail below for Ra = 0. 



(52) 






(53) 



We note that the dependence on the Bond number 
Bo enters through the quantity [p*-l]Bo-a^. 
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For the values given in Table I, the additional constants 
are given by 

C3 =7.0898, c, =4.6055x10"', c, =2.6859. (55) 

Note that if Cr = in Eq. (51) the two-term approxima- 
tion, 



Ma=^a-'+^, 



(56) 



is independent of Bo. This is a result of the boundary 
condition (29); for Cr = 0, the interface is rigid and the 
Bond number has no effect on the system. 

7.2 Large Wavenumbers 

We next consider the limit of large wavenumbers for 
a stationary mode. For our system the numerical results 
suggest that buoyancy effects and interface deforma- 
tion are unimportant in this limit, so we also consider 
the formal limit of small Crispation number, Cr -^ 0, 
along with Ra = 0. For Ra = 0, the governing equations 
for the velocity field are decoupled from the thermal 
field, which simplifies the analysis. For Cr — > the 
dimensionless form of the normal momentum balance, 

Cr(^"-^^)-Bo(p*-l)^ + fl'^=2Cr(^*w;-ft)^^, (57) 

then reduces to 

Bo(p*-l)^-fl'^ = 0. (58) 

For Bo (p* - 1) - a^ ;^ 0, we conclude that the interface 
deformation vanishes, ^ = 0. 

In the limit of large wavenumber the disturbances are 
concentrated near the interface and the effects of the 
outer boundaries are insignificant. The appropriate 
solution can then be computed by applying decay con- 
ditions in an unbounded domain as z -> ± 00. The verti- 
cal components of the velocity field are given by 

(59) 
The temperature fields are given by 



2a K 4a K 



G^'B 



(60) 



AaK 



« r 2 2 n -az , >-^ -1 

— [a z ]e +C^e 



la' 



4a' 



rp / \ Pr -I az Pr n«z, P r 2 2-i az , y^ az 

(z)= — -[az]e -[az]e -^^[a z ]e +C„e . 



Ad 



The corresponding horizontal velocities and pressures 
are 



u''=i{B^-AJe--iB^[az-\e-'\ 
u^ =i{B^+A^)e''' +iB^iazY\ 



and 



The interfacial boundary conditions are 
of =0, Q)^ =0, w" =w^ 

fa ^f^^ j^*fa_fp ^^^ 



(62) 



(63) 



(64) 



(65) 



(AiX" -u^,) + ia{lxocr -d)^)-/aMar" =0. (66) 

The interface boundary conditions determine the re- 
maining six constants A^, B^, C^^^A^^B^, and C^. 

The velocity boundary conditions require A^ = A^ = 0, 
and B^ = Bp, and the solution to the thermal problem 
then gives 



C=C, 



(1-K-*) 



4a'K-*(l + F) 



B^ 



(67) 



The remaining boundary condition gives the dispersion 
relation, 



Ma=8a'(l+/i )k-' 



., .(1 + r) 



O-K') 



(68) 



(61) 



For the values in Table I this gives Ma = 44.2 a^. 
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